In my experience, get_jacobian is slow because forming the Pynumero NLP is slow. If we have flowsheet-level variables, we end up getting a Jacobian for the flowsheet and every subcomponent. I'm not sure how slow forming the NLP is going to be for PropertyBlock 43 in the MEA column, but you're still at least doubling the work you're doing.

I would suggest getting the Jacobian and NLP once at the beginning of constraints_by_jacobian_norm, then just referring to rows and columns of the overall object as you want to scale constraints.

I have been thinking about there best way to do this, and did realize I was probably getting the Jacobian more times than necessary. I was not aware of the issue of having higher-level components in the problem causing the entire problem to be emitted (and given how properties work that will be an issue). I am also thinking there might need to be some checks to deal with deactivated components (and that it might be possible to deactivate unnecessary components before getting the Jacobian to reduce the problem size).

I think I've fixed this so this method will only be called once for any call to the scaling method.

constraints_by_jacobian_norm (line 156) collects all the constraints we want to scale and the top level block, and then only calls the autoscaler once.

Based on your comment, should we consider adding another method for variable-based norm scaling?

I do not know if this makes sense or not - comments would be welcome (@dallan-keylogic ).

If we scale variables first based on order-of-magnitude, we can then scale constraints based on row/column norm. However, norm-based variable scaling probably won't be useful for a fully unscaled model, because the most common Jacobian entries will be O(1).

For a partially-scaled model, though, filling in variable scaling factors based on column norms might be useful.

It'd be nice to have a method like report_scaling that would give you vars and constraints in one column, scaling factor values in the 2nd column, and the scaled value in the third column.

The first two you can get by doing suffix.display(). The third part - having scaled value is easy enough for Vars, but what is the "scaled value" for a Constraint? Scaled residual is easy to get, but not always that meaningful (as the scaled Jacobian is often more meaningful).

Actually, adding to my other suggestion on thinking about a reporting method for Autoscaler, these would also be good to show in report output (e.g., condition number, extreme jac rows/columns or at least the number of each).

I realize that the DiagnosticsToolbox already has the capability to report on some of these things, but just thinking it'd be nice to be able to apply scaling via AutoScaler and subsequently check what you really did to model scaling with a report method.

That is what the Diagnostics tools are for - putting those methods here would just be duplication of code (and the diagnostics also tell you a lot more). The eventual documentation will highlight that you really need to use diagnostics, scaling and initialization tools together to get the best results.

Will the diagnostics toolbox reference these new tools?

@agarciadiego I am not sure there - we probably need to think about how and where we draw the line between the two. Up until now, the Diagnostics Toolbox has been about finding the issues, and hasn't said much about how to fix them (as that is not so easy to explain concisely). We can definitely mention it in the docs, and maybe we could have the display methods mention tools to help fix issues (where appropriate).

These are all issues conditionally, although extreme Jacobian entries is more conditional than the others.
I believe "constraints with potential cancellations" aren't typically scaling issues but model formulation issues. "variables with values close to zero" is redundant with "variables with very large or very small values". "Variables with values close to bounds" typically end up being variables close to zero. I don't recall an example with a variable close to a nonzero bound that turned out to be bad scaling.

You could rank order them:

1. Variables and constraints with extreme Jacobian column/row norms
2. Variables with very large or very small values
3. Gap
4. Variables close to bounds
5. Constraints with mismatched terms
6. Extreme Jacobian Entries

As a side note, an instructive exercise for constraints with mismatched terms is to imagine setting the mismatched term to zero. If the Jacobian is nonsingular, it's probably fine. If it becomes singular, it's probably an issue.

For this documentation I deliberately didn't go into too much detail. There is definitely room for a more detailed discussion of all these issues, but for the basic user I just wanted to give a list of things to look for.

We could use the actual submodel objects if we did a ComponentMap. I assume you're aware of the option and elected for strings and dictionaries for some performance reason.

I think the reason for this was that I was thinking we did not have the objects yet, but as this is a user-facing function (and not a set of defaults) that is not the case. I will update this too.

I was going to ask about what happens if multiple scaling suffixes for the same variable is encountered on different blocks, but it looks like the full block structure is recreated in this function.

For all of these tools, we assume that the IDAES Controlled Suffix will be on the parent block of the component. Pyomo does support putting scaling factors elsewhere, but that would require the user to manually set those (and thus they are responsible for managing them).

Is there a way to suppress debug messages like this for the scaling routines only? When I turn on Debug logging, it's typically to look at the IPOPT output of an initialization solve. I guess that involves passing the output level to the initialization object, so I imagine that the scaling routines take an outlevel argument in the same way.

Not quite - you should be able to control the logging level by each subpackage in IDAES. I.e. you can set different levels for idaes and idaes.core.scaling.

Are ConcreteModels instances of Blocks?

Yes.

This method is technically correct, but I'm not sure how useful this concept of a "nominal value" is going to be. We'll see how it gets used.

This is the underlying function that is used by the expression walker to get the nominal value of an expression/constraint. It is primarily intended to be used in the expression walker (or limited similar circumstances).

Once this gets merged, I want to see what's causing the condition number to still be so badly behaved.

My guess is the fact that the model is uninitialized, but I gave it good initial guesses for the actual final state - i.e. there is a mismatch between the current uninitialized state and the scaling factors it was given. If we went and solved the model, I suspect the value would be much lower.

Not for this PR, but schemes could be explained more in the documentation